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ABSTRACT 

We used cosmological simulation with self-consistent radiative transfer to investigate 
the physical nature of Lyman limit systems at z = 4. In agreement with previous studies, 
we find that most of Lyman-limit systems are ionized by the cosmological background, 
while higher column density systems seem to be illuminated by the local sources of 
radiation. In addition, we find that most of Lyman limit systems in our simulations 
are located within the virial radii of galaxies with a wide range of masses, and are 
physically associated with them ("bits and pieces" of galaxy formation). While the 
finite resolution of our simulations cannot exclude an existence of a second population 
of self-shielded, neutral gas clouds located in low mass dark matter halos ( "minihalos" ) , 
our simulations are not consistent with "minihalos" dominating the total abundance of 
Lyman limit systems. 

Subject headings: cosmology: theory - cosmology: large-scale structure - cosmology: Lyman- 
limit systems 



1. Introduction 

The hydrogen absorption system observed 
in spectra of distant quasars are tradition- 
ally subdivided into Lyman-a forest (iVm < 
1.6 x 10 17 cm 2 ), Lyman-limit systems (1.6 x 
10 17 cm 2 < Nm ^ 10 20 cm 2 ), and damped 
Lyman-a systems (iVm > 10 20 cm 2 ). Amaz- 
ingly, among these three classes, the Lyman- 
limit systems are least understood, both ob- 
servationally and theoretically. This is quite 
unfortunate, given that the Lyman-limit sys- 
tems dominate the absorption of ionizing pho- 
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tons in the universe (Miralda-Escude 2003), 
and are, therefore, crucial for understanding 
the transfer of ionizing radiation in the IGM 
and the interactions between the early galax- 
ies and quasars and their environments. 

Previous attempts to model Lyman-limit 
systems in cosmological simulations identified 
them with low mass dark matter halos, filled 
with neutral, self-shielded gas (Katz et al. 
1996; Gardner et al. 1997, 2001). On the 
other hand, the observations of neutral hydro- 
gen in the vicinity of the Milky Way galaxy 
uncovered a large population of High Veloc- 
ity Clouds (HVCs) that could plausible be 
low redshift counterparts of Lyman-limit sys- 
tems (Oort 1966; Verschuur 1969; Blitz et al. 
1999; Maloney k Putman 2003; Putman et al. 
2003). 

In this paper we use cosmological simula- 
tions with the self-consistent treatment of the 



transfer of ionizing radiation to address the 
question of the nature of Lyman-limit sys- 
tems, as their appear in the simulations. An- 
other advantage of the simulations we use here 
is that they are designed to model the reion- 
ization of the universe, and so achieve a rea- 
sonable agreement with the SDSS observa- 
tional data on the Lyman-a forest at z > 5 
(Fan et al. 2006; Gnedin & Fan 2006). 

2. Simulations 

Simulations used in this paper have been 
run with the "Softened Lagrangian Hydro- 
dynamics" (SLH) code (Gnedin 2000, 2004). 
Simulations include dark matter, gas, star for- 
mation, chemistry, and ionization balance in 
the primordial plasma, and three-dimensional 
radiative transfer. The simulations described 
in this paper use a flat ACDM cosmology with 
values of cosmological parameters as deter- 
mined by the first year WMAP data (Spergel 
et al. 2003). 1 

We use two sets of simulations in this pa- 
per. The two sets differ by the size of the 
computational box: in the first set the com- 
putational box has the size of 4/i _1 comoving 
Mpc, while in the second set the box is 8/i _1 
comoving Mpc on a side. 

The primary simulations in each of two sets 
are runs L4N128 and L8N128 from Gnedin 
& Fan (2006). Both of them include 128 3 
dark matter particles and the same number of 
quasi-Lagrangian mesh cells for the gas evo- 
lution, as well as a smaller number of stellar 
particles that formed continuously during the 
simulation. 

As Gnedin & Fan (2006) emphasized, these 
simulations provide an adequate fit to the evo- 
lution of the mean transmitted flux in the 
Lyman-a forest between 5 < z < 6.2, but 
overproduce the opacity in the forest at z < 5. 

1 Q. M = 0.27, Q.a = 0.73, h = 0.71, Q B = 0.04, n s = 1. 



As we show below, they also overproduce 
Lyman-limit systems. 

The simulations, however, contain ad- 
justable parameters; the one of importance 
here is the ionizing efficiency parameter euv, 
that controls the amount of ionizing radiation 
emitted per unit mass of stars formed in the 
simulation. This parameter depends on the 
numerical resolution of the simulation and on 
the details of stellar feedback and radiative 
transfer implementations, and has no clear 
physical interpretation. It can be thought of 
as a fraction of ionizing radiation escaping 
from the spatial scales unresolved in the sim- 
ulation (it has no direct relation to the widely 
used "escape fraction" quantity). 

Since we treat the ionizing efficiency as a 
free parameter, we can adjust it to obtain a 
better agreement with the observational data. 
For that purpose, we run additional simula- 
tions in which we increased the ionizing effi- 
ciency by an additional factor q. Because in 
this paper we only concentrate on comparing 
the simulations with the data at z = 4, and 
in order to save computational resources, we 
only rerun the simulations with different val- 
ues of the ionizing intensity from z = 4.3 to 
2 = 4. The redshift interval of Az = 0.3 corre- 
sponds to the time interval of 130 Myr, which 
is some 4,000 larger than the mean photoion- 
ization time at this redshift (McDonald et al. 
2002), and, thus, should be more than suf- 
ficient for establishing a new photoionization 
equilibrium even in the lowest density regions. 
To be on the safe side, we also rerun one of 
the simulations from z = 4.6 to z = 4, and 
the results for the Lyman-a mean transmit- 
ted flux and the column density distribution 
of the Lyman-limit systems from that simula- 
tion are indistinguishable from the analogous 
run started at z = 4.3. 

In order to make references to specific sim- 
ulations transparent, we label each simula- 
tion with a letter L followed by the value of 
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the linear size of the computational volume 
(measured in h~ l Mpc), followed by the letter 
"q" and the value of the ionization enhance- 
ment parameter q. For example, the primary 
simulations in each set are thus labeled L4ql 
and L8ql, while a 4/i _1 Mpc simulation with 
6 times higher ionizing efficiency is labeled 
L4q6. In the 4/i _1 Mpc set we run simula- 
tions L4ql, L4q2, L4q6, and L4ql0, while in 
the 8/i _1 Mpc set we run L8ql and L8q6. 

The Lyman Limit systems are found by 
casting 10000 lines of sight through the sim- 
ulation box, allowing us to sample the sim- 
ulation box finely enough to find even rare 
absorbers. The lines of sight are cast in ran- 
dom directions and cover a redshift range of 
5z = 0.033 for the 8/1" 1 Mpc boxes. This nar- 
row extent in redshift space makes it difficult 
to search for the Lyman-limit edge in the spec- 
tra. Thus, to find the regions of high optical 
depth along the lines of sight, we search for 
peaks in Xdr^jdX. 

This is similar to searching for an absorp- 
tion line. For example, a Lyman-a absorption 
line profile can be described with (we ignore 
natural width of the line here): 




When considering Xdrn/dX instead of tll, 
the Lyman limit edge becomes the delta- 
function in frequency, in full analogy with the 
absorption line profile. Thus, a Lyman-limit 
system appears as a peak in the quantity 

1 f , ,(^') 2 ,dT LL 
TLL = b^/n J ~~d\ r ' (2) 

This substitution allows us to use soft- 
ware tools, already developed for finding ab- 
sorption lines in synthetic spectra, for find- 
ing Lyman-limit systems in the simulated syn- 
thetic spectra. In particular, the value of 
tll at the "line" center is directly propor- 
tional to the neutral hydrogen column density, 



tll = ctllNhi, where a L L = 6.3 x 10 18 cm 2 
is the hydrogen photoionization cross-section. 
All regions in the spectra where tll > 1-0, 
corresponding to column densities of Nhi > 
1.6 x 10 17 cm -2 , are selected as Lyman-limit 
systems. 

After casting the lines of sight, the posi- 
tions of these regions (in 3D space) and their 
opacities are determined. We also compute 
other characteristics such as local photoion- 
ization rate, column density and neutral frac- 
tion of the gas. This information allows us 
to investigate whether the systems are ion- 
ized by local sources or by the cosmological 
background radiation. 

Finding the positions of the Lyman-limit 
systems in the simulation box and determin- 
ing the position of the galaxies in the sim- 
ulation allows us to determine whether each 
Lyman Limit system can be associated with a 
galaxy. 

Galaxies in the simulation are identified 
with gravitationally bound objects, found 
with the DENMAX algorithm (Bertschinger 
& Gelb 1991). The DENMAX algorithm 
works by constructing the finite resolution 
smooth total (dark matter + gas + stars) 
density field from the simulation data. It then 
identifies density maxima and associates them 
with bound objects. This approach has a ma- 
jor advantage for our purposes here, since it 
allows us to assign a specific meaning to the 
word "associated" . Since each galaxy is repre- 
sented by DENMAX as a density hill, we will 
call all points in space that lie on that hill 
as "associated" with that galaxy. In other 
words, a point of space is "associated" with a 
given galaxy if, by moving against the density 
gradient from that point one would eventually 
end up at the center of the galaxy. 

DENMAX separates all points of space into 
four distinct categories: 

• points gravitationally bound to a given 
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galaxy; 



• points associated with a given galaxy 
but not gravitationally bound to it (i.e. 
points lying on the galaxy density hill, 
but outside the virial radius); 

• unassociated points (i.e. points not as- 
sociated with any galaxy, like points in 
the middle of the void); and 

• points associated with unresolved galax- 
ies (i.e. points belonging to the density 
hill too small to be qualified as a re- 
solved galaxy). 

We use this classification below to describe 
locations of Lyman-limit systems in space. 

3. Results 

The main goal of this work is to investigate 
the physical properties of the Lyman limit sys- 
tems themselves, such as neutral fraction and 
photoionization rate, but also what kind of 
galaxies they are associated with. 

To begin, we have to normalize the amount 
of photons in the simulation box to the ob- 
servations, so that the population of simu- 
lated Lyman-limit systems corresponds to the 
systems that are observed. This is done us- 
ing the number density of Lyman limit sys- 
tems per column density and adjusting the 
ionization efficiency to reduce the amount of 
Lyman-limit systems produced. Another ob- 
servational parameter that can be used to fit 
the simulations to observations is the mean 
transmitted flux in Lyman-a. 

Figure 1 shows the number density of the 
Lyman limit systems versus their column den- 
sity for 3 different runs with varying ionization 
efficiency q. The top run, with its fitted line 
overlaid, corresponds to L4ql, which is the the 
run with the default value of the ionization ef- 
ficiency. The middle run is L4q6, where the 
ionization efficiency is increased by a factor of 
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Fig. 1. — Shown are the number densities of 
Lyman-Limit systems for three different runs, 
each with 4/j -1 Mpc box size. The fits are 
overlaid. Top: L4ql, middle: L4q6, bottom: 
L4ql0. 

6 and the bottom run is run L4ql0, with a 
ionization efficiency q = 10. The overlaid fits 
are best-fit power laws with the slope —1.5 
(Petitjean et al. 1993; Miralda-Escude 2003). 

For the L4q6 run the fitted line exactly cor- 
responds to the observed number density ob- 
tained by Storrie-Lombardi et al. (1994). The 
other two fits have the same power law, but 
different amplitudes. L4ql had a number den- 
sity of Lyman-limit systems too high to fit the 
observations. This illustrates that the original 
simulations, while fitting the mean transmit- 
ted flux in the Lyman-a forest at z > 5, sig- 
nificantly overpredict the amount of Lyman- 
limit systems at z = 4. 

Figure 2 similarly shows the number den- 
sity distribution with column density for the 
two 8/i -1 Mpc runs L8ql and L8q6 and for 
comparison the 4/i _1 Mpc run L4q6. The 
highest distribution is again L8ql with the 
lowest ionizing efficiency. The other 8/i -1 Mpc 
run also has the same q = 6. There is good 
agreement between the 4 and the Sh~ l MpC 
run at Nhi > 10 17 cm 2 . At lower column den- 
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Fig. 2. — Similarly to Figure 1, this graph 
shows the number density distribution of Ly- 
man limit systems for two 8/i _1 Mpc runs (in 
black, L8ql on top and L8q6 on bottom) and 
4/i _1 Mpc run L4q6 (in grey). Their power 
law fits are included. 

sities the larger box size lacks enough resolu- 
tion to fully resolve the Lyman-a forest. Thus 
the number of Lyman-a forest systems is too 
low at lower column densities, which can be 
seen from the disagreement for the two q = 6 
distributions. 

In Figure 3 the dependence on the param- 
eter q is illustrated. The solid lines show the 
mean transmitted flux in the Lyman-a for- 
est normalized by the observational values of 
0.55 at z = 4. For higher q there is more flux 
transmitted as the increased ionizing intensity 
ionizes more of the initially neutral gas. The 
dashed lines show the amplitude of the power 
law that is needed to fit the number density 
of Lyman-limit systems as described above in 
Figures 1 and 2. Here, for lower ionization 
efficiency the offset is lower, for q = 6 the off- 
set is unity, since the number density distribu- 
tion is best fit by the observational value. The 
horizontal line shows the observational values 
for both parameters since we normalize by the 
observed value. The gray lines correspond to 
a 20% error, as is expected from the observa- 
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Fig. 3. — This graph shows the ionization effi- 
ciency q that was used as a fitting parameter 
for the simulations on the x-axis and shown 
on the y-axis are the two observational con- 
straints: amplitude of the number density fit, 
A, and the mean transmitted flux F mean . 

tions. This plot illustrates that we can fit two 
different observables with just one parameter, 
the ionization efficiency, and that increasing 
this parameter by a factor of 6 from the values 
used in Gnedin & Fan (2006) fits the observa- 
tions at z = 4. 

Figure 4 shows the distribution of optical 
depth of the Lyman-limit systems versus hy- 
drogen neutral fraction. There is a tight rela- 
tion of the optical depth with neutral fraction, 
the higher tll, the less ionized the system. 
The remarkable conclusion from this is that 
the systems are substantially ionized and are 
not neutral. This implies that the Lyman- 
limit systems detected in this simulation are 
not small self-shielded clumps of neutral gas, 
but rather diffuse highly non-spherical clouds 
of gas that happen to be oriented in their 
longest direction toward an observer. 

In Figure 5 it is similarly shown how the 
optical depth depends on the photoionization 
rate for each Lyman-limit system. These two 
variables are not as tightly correlated as the 
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Fig. 4. — Distribution of Lyman limit system 
optical depth versus the neutral fraction of 
gas. 



Fig. 5. — Similar to Figure 4, this graph shows 
the distribution of optical depth versus the 
normalized photoionization rate. 



optical depth and the neutral fraction. For 
lower tll systems the normalized photoion- 
ization rate T/(T), where T is the average T, 
is about unity. This means that these systems 
are ionized by the cosmological background 
and not by local sources. Systems with higher 
tll appear to be ionized by local sources be- 
cause their normalized photoionization rate is 
higher. This result is in agreement with pre- 
vious findings by Miralda-Escude (2005)and 
Schaye (2006). 

One of the questions we investigate is 
whether the Lyman-limit systems are asso- 
ciated with high or low mass galaxies. At 
first, we have tried to use a method similar 
to one used by observational studies of the 
environment of the Lyman-limit systems at 
low redshifts (Penton et al. 2002). In obser- 
vations, the Lyman-limit systems are often 
associated with a galaxy by finding the near- 
est object, either along a line of sight or near 
it. A similar method was also used to analyze 
the simulations in Katz et al. (1996). There, 
a relation between the distance to the near- 
est galaxy and Nm of the Lyman-limit sys- 
tems was shown, and it was concluded that 
the column density correlates inversely with 



projected distance, and that Lyman-limit sys- 
tems either lie in the outer parts of massive 
protogalaxies or closer to the center of less 
massive galaxies. 

In our case we found that associating the 
Lyman-limit systems with the nearest ob- 
ject is not necessarily the best option, be- 
cause such association is strongly dependent 
on the galaxy sample. When we decrease 
the mass limit of the sample, thus includ- 
ing less and less massive galaxies in the sam- 
ple, we find progressively lower mass galax- 
ies closer and closer to a typical Lyman-limit 
system. This, of course, does not imply that 
the Lyman-limit systems are physically asso- 
ciated with these low mass galaxies, because 
low mass galaxies are more numerous, and 
the distances between even a random subset 
of points in space to galaxies of progressively 
smaller masses will be progressively smaller. 

Thus, associating the Lyman-limit systems 
with nearest galaxies could lead to the con- 
clusion that most of the Lyman-limit systems 
are located in small galaxies. However, us- 
ing the physical association as permitted by 
the DENMAX algorithm (§2) - i.e. associating 
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a Lyman-limit system with the galaxy whose 
"density hill" (i.e. dark matter halo) the 
Lyman-limit is located on - we find a different 
result: the Lyman-limit systems are associ- 
ated with galaxies in a wide range of masses 
and the small galaxies near them are often 
physically unrelated satellites of the larger 
galaxy. 
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Fig. 6. — Distribution of optical depths tll 
of the Lyman-limit systems versus the mass 
of their associated galaxy. 

Figure 6 shows the masses of the galaxies 
that the Lyman-limit systems are associated 
with versus the systems' optical depth. There 
is no strong dependence on galaxy mass for 
the optical depth, meaning that the size of the 
galaxy does not appear to strongly influence 
the opacity of the Lyman-limit system. 

The bottom panel of Figure 7 is a graph 
of the distance of the Lyman-limit system to 
its associated galaxy and the mass of that 
galaxy. Overlaid on this correlation is the dis- 
tance that corresponds to the virial radius for 
each galaxy mass. The resolution limit of our 
simulation is about 0.001 Mpc, so that we re- 
solve all the distances to the galaxies. The 
top panel shows the number of galaxies per 
mass bin for the same mass range as the bot- 
tom panel. It is clear from this histogram that 
there are many more small galaxies that have 
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Fig. 7. — Top Panel: Number of galaxies per 
mass bin. Bottom Panel: Distribution of dis- 
tances of the Lyman-limit systems from their 
associated galaxy with respect to the mass of 
that galaxy. 

no Lyman-limit associated with them. 

Combining the information from this graph, 
we find that most of the Lyman-limit systems 
are within the virial radius of a galaxy, often a 
relatively massive one. This means that they 
are either in orbit around or infalling into a 
galaxy, and not isolated clumps of neutral gas. 
Using the DENMAX classification of points in 
space discussed in §2, we find no Lyman-limit 
systems that are unassociated or associated 
with unresolved galaxies. 

Figure 8 extends the discussion from the 
previous graph by showing the number of 
Lyman-limit systems associated with galaxies 
of different masses. There are many more sys- 
tems per galaxy in the high galaxy mass range 
around M gal ~ 10 10 to M gai ~ 10 11 than for 
lower masses. Since dN ga i/d log(M) ~ M _1 , 
the number of galaxies sharply increases at 
low masses. In Figure 8 the overlaid power law 
(oc M 1 ) shows the equal number of Lyman- 
limit systems per logarithmic bin in mass of 
galaxies they are associated with. The figure 
shows that only a small fraction of all Lyman- 
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Fig. 8. — Number of Lyman-limit systems per 
associated galaxy. Overlaid is the fit corre- 
sponding to an equal amount of Lyman-limit 
systems per unit log in galaxy mass. 

limit systems are associated with galaxies of 
M ga i < 10 10 M Q . However, this effect might 
also be related to resolution effects, meaning 
that the Lyman-limit systems are associated 
with galaxies of all mass ranges equally. 

4. Conclusions 

We show that the column density distribu- 
tion of Lyman-limit systems at z = 4 in the 
cosmological simulations that agree with the 
SDSS measurements of the Lyman-a absorp- 
tion at z > 5 in the spectra of high redshift 
quasars has the same shape as the observed 
distribution. The abundance of the Lyman- 
limit system is not necessarily in agreement 
with the data, but we manage to achieve 
the agreement with both the observed column 
density distribution of the Lyman-limit sys- 
tems and with the mean transmitted flux in 
the Lyman-a forest by adjusting a single pa- 
rameter - ionizing efficiency - at z < 5. 

We find that, in the simulation that 
agrees with the observational data, the low- 
est column density Lyman-limit systems are 
mainly illuminated by the cosmological back- 



ground, in agreement with previous findings 
of Miralda-Escude (2005) and Schaye (2006). 

However, we also find that all Lyman-limit 
systems that are resolved in our simulations 
(Njn < 10 20 cm 2 ) are highly ionized and are 
physically associated with (are in orbit around 
or falling into) the gravitational well of a large 
range of galaxy masses (M tot > 1O 1O M ). In 
other words, they appear as residual "bits and 
pieces" of the galaxy formation process rather 
than self-shielded, highly neutral dark matter 
"minihalos" (M tot < 1O 9 M ). 

Because of the finite spatial and mass res- 
olution of our simulations, we can not ex- 
clude an existence of the second population 
of Lyman-limit systems composed of "miniha- 
los" . However, the fact that we are able to fit 
both the Lyman-limit column density distri- 
bution and the mean transmitted flux in the 
Lyman-a forest at z = 4 with one adjustable 
parameter may indicate that the population of 
Lyman-limit systems arising in "minihalos" is 
not the dominant one. Other studies have also 
shown that we can reproduce the observed 
Lyman-a forest, so that the hidden popula- 
tion cannot have any influence on the Lyman- 
a forest. Cross-correlating Lyman-limit sys- 
tems and galaxies could also help determine 
whether this population is related to "mini- 
halos". 
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